In a city as diverse as Boston, different communities may be less prepared for the extreme effects of climate change. Factors such as older or younger age, pre-existing medical conditions and low-income can place a person at greater risk to extreme heat.

With the urban heat island effect, urban areas, lacking in greenery, are more susceptible to extreme heat than rural ones. Extreme heat can accumulate on heat-absorbent surface areas, such as buildings and pavement. But within Boston itself, neighborhoods vary by the amount of open space and trees, which contributes to varying extents of urban heat in different communities on the same hot summer day, as seen in Metropolitan Area Planning Council’s map for the Boston metropolitan area. .

This year, Boston had its hottest summer on record. As the Earth’s climate continues to warm, identifying which neighborhoods are more at risk to extreme heat is crucial for climate preparedness.

This data project, told in three parts and developed with R, will help provide visualizations of the disproportionate effect of a warmer Boston for specific communities and examine how historical redlining may have had a lasting effect on those communities today.

Source: Climate Ready Boston

Source: Mapping Inequality & Analyze Boston A preview of visualizations to come

Setup

# Load required libraries
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.5     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.0.2     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(htmlwidgets)
library(forcats)

Loading data

## Rows: 180 Columns: 16
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): Name
## dbl (15): FID, GEOID10, AREA_SQFT, AREA_ACRES, POP100_RE, HU100_RE, TotDis, ...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Preparing Data

# use dplyr to rename and select certain columns
climate_ready_sel <- climate_ready %>% select(FID, GEOID10, tract_population = POP100_RE, tract_housing = HU100_RE, neighborhood = Name, low_income = Low_to_No, LEP, older = OlderAdult, POC = POC2, med_illness = MedIllnes, children = TotChild) %>% 
  glimpse()
## Rows: 180
## Columns: 11
## $ FID              <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16…
## $ GEOID10          <dbl> 25025010405, 25025010404, 25025010801, 25025010702, 2…
## $ tract_population <dbl> 5522, 5817, 2783, 2400, 3173, 3059, 4804, 2791, 4985,…
## $ tract_housing    <dbl> 994, 1862, 1899, 1643, 1283, 2106, 2136, 1418, 2726, …
## $ neighborhood     <chr> "Mission Hill", "Fenway", "Back Bay", "Back Bay", "Fe…
## $ low_income       <dbl> 1191, 2387, 72, 187, 895, 309, 1024, 1740, 1967, 1560…
## $ LEP              <dbl> 1522, 2443, 462, 472, 931, 737, 1406, 1771, 2804, 167…
## $ older            <dbl> 331, 56, 390, 285, 36, 428, 382, 31, 837, 113, 16, 96…
## $ POC              <dbl> 1755, 1749, 447, 320, 1039, 517, 1664, 1083, 1123, 19…
## $ med_illness      <dbl> 2131.22, 2201.14, 1214.76, 1014.20, 1181.78, 1245.52,…
## $ children         <dbl> 60, 77, 281, 86, 13, 232, 85, 31, 35, 142, 34, 168, 0…
# sum different population data columns based on neighborhood
climate_ready_analysis <- climate_ready_sel %>% group_by(neighborhood) %>% 
  summarize(neighborhood_population = sum(tract_population), POC = sum(POC), LEP =sum(LEP),
            older = sum(older), med_illness = sum(med_illness), low_income = sum(low_income), children = sum(children)) %>% glimpse()
## Rows: 23
## Columns: 8
## $ neighborhood            <chr> "Allston", "Back Bay", "Bay Village", "Brighto…
## $ neighborhood_population <dbl> 7592, 28634, 10850, 60821, 16439, 69695, 40517…
## $ POC                     <dbl> 3011, 6923, 6008, 19768, 3981, 50778, 25459, 9…
## $ LEP                     <dbl> 3444, 9075, 5640, 22914, 5968, 29284, 17845, 7…
## $ older                   <dbl> 144, 3759, 1711, 5847, 1811, 6535, 4147, 496, …
## $ med_illness             <dbl> 2831.85, 11842.96, 4292.93, 23936.80, 6461.47,…
## $ low_income              <dbl> 3300, 5316, 3929, 17067, 4157, 22749, 13698, 6…
## $ children                <dbl> 173, 1686, 1249, 4244, 3301, 17424, 8665, 453,…
view(climate_ready_analysis)

#write.csv(climate_ready_analysis, "data/climate_ready_analysis.csv")

Checking Boston neighborhoods

# Checking neighborhood names in Boston Neighborhoods dataset

# Load CSV for Boston Neighborhoods
boston_neighborhoods <- read_csv("data/Boston_Neighborhoods.csv")
## Rows: 26 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Name
## dbl (6): OBJECTID, Acres, Neighborhood_ID, SqMiles, ShapeSTArea, ShapeSTLength
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
boston_neighborhoods %>% glimpse() # count: 26
## Rows: 26
## Columns: 7
## $ OBJECTID        <dbl> 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40…
## $ Name            <chr> "Roslindale", "Jamaica Plain", "Mission Hill", "Longwo…
## $ Acres           <dbl> 1605.56824, 2519.24539, 350.85356, 188.61195, 26.53984…
## $ Neighborhood_ID <dbl> 15, 11, 13, 28, 33, 27, 26, 14, 16, 32, 2, 8, 4, 31, 3…
## $ SqMiles         <dbl> 2.51, 3.94, 0.55, 0.29, 0.04, 0.02, 0.12, 0.20, 3.29, …
## $ ShapeSTArea     <dbl> 69938272.9, 109737890.8, 15283120.0, 8215903.5, 115607…
## $ ShapeSTLength   <dbl> 53563.913, 56349.937, 17918.724, 11908.757, 4650.635, …
summary_neighborhoods <- boston_neighborhoods %>% select(Name) %>% 
  arrange(Name) %>% glimpse()
## Rows: 26
## Columns: 1
## $ Name <chr> "Allston", "Back Bay", "Bay Village", "Beacon Hill", "Brighton", …
# Noting differences between neighborhood names
view(summary_neighborhoods) 
# Additional neighborhoods: Beacon Hill, Chinatown and Downtown
# Longwood

view(climate_ready_analysis$neighborhood) # count: 23
# Longwood Medical Area

# download csv of neighborhood names
# write.csv(summary_neighborhoods, "data/neighborhoods.csv")

Part 1: A Heat Map of Population Distribution of Vulnerable Communities

In order to examine different populations in Boston neighborhoods, I used a 2016 dataset of vulnerable social groups from Climate Ready Boston, the city’s initiative for preparing for the long-term effects of climate change. This dataset was obtained from Analyze Boston.

The original dataset, derived from research by certified emergency manager Atyia Martin, has population data organized by U.S. Census tract for older adults, children, people with limited English proficiency (LEP), low-income, people of color and individuals with medical illnesses, such as asthma, heart disease, emphysema, bronchitis, cancer, diabetes, kidney disease and liver disease. It also includes the neighborhood name and total populations for each Census tract.

Some things to note when interpreting this data: * The population data is based on 5-year estimates from the U.S. Census American Community Survey from 2008-2012, which may be somewhat outdated and not directly reflect population data this year. * Many of these same individuals who are part of one group may be counted toward another or more groups, but the data does not identify population of people in multiple groups.

Making a heat map

In order to simplify the data for visualizing, I used dplyr. I found the totals for each group across neighborhood names. Note that the original dataset only has 23 neighborhoods and does not include Beacon Hill, Chinatown or Downtown. Knowing I wanted to create a heat map, I converted the data into long-form format using the tidy function in the tidyr package.

Then, using ggplot2, I created this heat map to visualize the percent makeup of each population group compared to the total neighborhood population. This allows for some standardized comparison across different neighborhoods and you can identify which communities may be in need of greater climate preparedness, depending on the frequency of deeper tiles.

The following interactive heat map was created with the plotly package and its ggplotly function. Hover over the tiles and a tooltip will pop up.

# Reformat the data from wide to long dataframe format using "tidyr" package (included in tidyverse). See "Reshape Data": https://github.com/rstudio/cheatsheets/blob/master/data-import.pdf.
# Add a percentage column
climate_ready_percents <- gather(climate_ready_analysis, "POC", "low_income", "LEP", "older", "children", "med_illness", key ="population_type", value= "population") %>% 
  mutate(population = round(population, digits = 0)) %>% 
  mutate(percentage = round((population/neighborhood_population), digits=2)) %>%
  glimpse()
## Rows: 138
## Columns: 5
## $ neighborhood            <chr> "Allston", "Back Bay", "Bay Village", "Brighto…
## $ neighborhood_population <dbl> 7592, 28634, 10850, 60821, 16439, 69695, 40517…
## $ population_type         <chr> "POC", "POC", "POC", "POC", "POC", "POC", "POC…
## $ population              <dbl> 3011, 6923, 6008, 19768, 3981, 50778, 25459, 9…
## $ percentage              <dbl> 0.40, 0.24, 0.55, 0.33, 0.24, 0.73, 0.63, 0.36…
# Heat Map, reference: https://www.r-graph-gallery.com/79-levelplot-with-ggplot2.html
# Add new column for plotly tooltip "text" in heatmap later
climate_ready_percents <- climate_ready_percents %>% 
  mutate(description = paste0("neighborhood: ", neighborhood, "\n","proportion: ", percentage, "\n", "population: ", population, "\n", "neighborhood population: ", neighborhood_population)) 

# Create heat map
# RColorBrewer color palettes: https://www.r-graph-gallery.com/38-rcolorbrewers-palettes.html
heatmap.climate_ready <- ggplot(climate_ready_percents, mapping = aes(x = population_type, y=neighborhood, fill = percentage, text = description)) +
  geom_tile(width=1, color = "gray") +
  scale_fill_distiller(palette= "YlOrRd", direction= +1) +
  labs(titles = "Boston Communities Vulnerable to Climate Change", subtitle = "Based on data from 2008-2012 American Community Survey 5-year Estimates, U.S. Census", x = "population", y = "neighborhood", fill = "proportion", caption = "Source: Climate Ready Boston") +
  theme_classic() +
  theme(axis.text.x = element_text(size = 8),
        axis.text.y = element_text(size = 8),
        plot.title = element_text(size=12, face = "bold", hjust = 0.5),
        plot.caption= element_text(size=7),
        plot.subtitle = element_text(size=9, hjust =0.5),
        legend.text = element_text(size = 6),
        panel.grid.major = element_line(size = 0.2, color = "gray"))
heatmap.climate_ready

#ggsave("heatmap.png", width=7, height = 5, unit="in")

## Make heatmap interactive using ggplotly
iheatmap.climate_ready <- heatmap.climate_ready %>% ggplotly(tooltip="description")
iheatmap.climate_ready
## Save html widget using library(htmlwidgets)
#saveWidget(iheatmap.climate_ready, "img/ggplotlyHeatmap.html")

Part 2: Neighborhood Comparisons of People of Color and Low-income Populations

While examining percentages of population is useful, it is equally important to identify which communities have the most vulnerable people by overall population. I chose the groups - people of color and low-income populations - to analyze more closely.

Making a Grouped Bar Chart

Using dplyr again, I filtered the same long-form dataset, derived from the Climate Ready Boston dataset, to only these two groups and reordered the neighborhoods by population values and then created a grouped bar chart using ggplot2.

The following grouped bar chart shows that Roxbury and Dorchester have the most people of color and low-income populations of the Boston neighborhoods. Note that these neighborhoods are also among the largest Boston neighborhoods.

# Reorder values displayed in graph using forcats library, see reference: https://www.r-graph-gallery.com/267-reorder-a-variable-in-ggplot2.html
# Grouped bar chart for low income and POC only - but by population
plot.POC_lowincome <- climate_ready_percents %>%
  filter(population_type %in% c("POC", "low_income")) %>%
  mutate(neighborhood = fct_reorder(neighborhood, population)) %>% ggplot(mapping = aes(x= population, y=neighborhood, fill= population_type)) +
  geom_bar(position = "dodge", stat = "identity", width =.75, alpha = 0.8) +
  scale_fill_manual(values = c("#ff8100", "#6821d8")) +
  labs(title = "Populations of Low-income and People of Color in Boston", x="population", y="neighborhood", fill = "population type", subtitle="Based on data from 2008-2012 American Community Survey 5-year Estimates, U.S. Census",
       caption = "Source: Climate Ready Boston") +
  theme_dark() +
  theme(axis.title = element_text(size=9), 
        axis.text = element_text(size = 7, color = "black"),
        axis.line = element_line(size = 1, color ="#747272"),
        axis.ticks.length.y = unit(.15, "cm"),
        axis.ticks.length.x = unit(.15, "cm"),
        plot.title = element_text(size=12, face="bold", hjust=0.5),
        plot.caption = element_text(size=7),
        plot.caption.position = "plot",
        plot.subtitle = element_text(size=9),
        legend.title = element_text(size = 9),
        legend.text = element_text(size = 7))
plot.POC_lowincome

#ggsave("img/POC_lowincome_barchart.png", width=8, height = 5, unit="in")

Part 3: How the Racist Past Affects Boston Today

Examining Racist Historical Housing Policy, Open Space and Low-income Populations

Researchers suggest that in the United States, people of color tend to be at greater risk to extreme heat and other public health emergencies, including the current COVID-19 pandemic. This is likely due to more people of color having susceptible health conditions and/or living in low-income houses.

However, researchers believe that historical redlining, a racist housing policy conducted in the 1930s by the federal agency Home Owners’ Loan Corporation, has lasting underlying effects, partly responsible for this racial disproportion, as reported by The New York Times. HOLC appraisers graded different areas from A (green), B (blue), C (yellow) and D (red) and in turn determined who would be allowed to receive federal and bank housing loans. The D-grade, a redlined area, would be the lowest grade and often assigned to neighborhoods with a majority of Black residents and other people of color.

Making choropleths

To visualize the effects of this for Boston, I created several choropleths using the ggplot2 function geom_polygon. I primarily worked with geoJSON files for this part of the project. To read each geoJSON in R, I used the geojsonio library and converted each dataset (fortified) using the broom package.

Historically redlined areas and present day Boston neighborhood borders

This map shows historically redlined areas and layers modern Boston neighborhood borders.

I obtained a geoJSON of 1938 HOLC map of Boston from Mapping Inequality, a project that has digitized historical redlining maps for major U.S. cities.

I also used a geoJSON dataset of Boston neighborhoods, obtained from Analyze Boston, in order to have neighborhood borders on the map. This dataset includes 26 neighborhoods and does include Beacon Hill, Chinatown or Downtown, as opposed to the earlier Climate Ready Boston dataset.

After plotting the redlined areas and layering the neighborhood borders, I found it quite difficult to label the neighborhoods with their respective names. As a workaround, I isolated a list of 26 neighborhoods and created a new dataset with the center of each neighborhood by longitude and latitude coordinates obtained from Google Search. This dataset is available in this Google Sheet. I then used geom_text_repel from the ggrepel package to plot these coordinates and neighborhood names as labels.

library(geojsonio)
## Registered S3 method overwritten by 'geojsonsf':
##   method        from   
##   print.geojson geojson
## 
## Attaching package: 'geojsonio'
## The following object is masked from 'package:base':
## 
##     pretty
library(broom)
library(ggrepel)
library(mapproj)
## Loading required package: maps
## 
## Attaching package: 'maps'
## The following object is masked from 'package:purrr':
## 
##     map
## Warning in RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid"): Ring Self-
## intersection at or near point -71.125747950000004 42.272347189999998
## SpP is invalid
## Warning in rgeos::gUnaryUnion(spgeom = SpP, id = IDs): Invalid objects found;
## consider using set_RGEOS_CheckValidity(2L)
## Warning in RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid"): Ring Self-
## intersection at or near point -71.10476122 42.272897690000001
## Warning in RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid"): Ring Self-
## intersection at or near point -71.132332660000003 42.293209140000002
## Warning in RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid"): Ring Self-
## intersection at or near point -71.066056529999997 42.373780060000001
## Warning in RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid"): Ring Self-
## intersection at or near point -71.126547799999997 42.371628090000002
## Warning in RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid"): Ring Self-
## intersection at or near point -71.077928909999997 42.34569887
## SpP is invalid
## Warning in rgeos::gUnaryUnion(spgeom = SpP, id = IDs): Invalid objects found;
## consider using set_RGEOS_CheckValidity(2L)
## Rows: 26 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Name
## dbl (3): Id, long, lat
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Part 3 - R analysis

# Choropleth in ggplot, reference: https://www.r-graph-gallery.com/327-chloropleth-map-from-geojson-with-ggplot2.html

# Choropleth 1: Historically redlined areas and present day Boston neighborhood borders
# Add new id2 column to boston_redlining_fort for simplified neighborhood grade (eg. A1 -> A) 
# Extract only first letter of id using base R substr function
# Source: Joachim Schork of Statistics Globe https://statisticsglobe.com/r-extract-first-or-last-n-characters-from-string
boston_redlining_fort <- boston_redlining_fort %>% 
  mutate(id2 = substr(id, 1, 1)) %>% glimpse()
## Rows: 1,639
## Columns: 8
## $ long  <dbl> -71.11985, -71.11990, -71.12110, -71.12184, -71.12292, -71.12452…
## $ lat   <dbl> 42.32301, 42.32178, 42.32161, 42.32102, 42.32026, 42.31958, 42.3…
## $ order <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 1…
## $ hole  <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, F…
## $ piece <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ group <fct> A1.1, A1.1, A1.1, A1.1, A1.1, A1.1, A1.1, A1.1, A1.1, A1.1, A1.1…
## $ id    <chr> "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1", "A1"…
## $ id2   <chr> "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A",…
#view(boston_redlining_fort)

## coordinates: rename "Name" column to "neighborhood"
coordinates <- data_coordinates %>% select(neighborhood = "Name", long, lat)
#view(coordinates) 
# This dataset has 26 neighborhoods, derived from Boston Neighborhoods csv

# Make Choropleth 1
# label regions using geom_text and boston_coord
# Use coord_map() for map-like rendering https://ggplot2.tidyverse.org/reference/coord_map.html
# Use geom_text_repel so neighborhood labels don't overlap on map. See: https://cran.r-project.org/web/packages/ggrepel/vignettes/ggrepel.html
plot.redlining <- ggplot() +
  geom_polygon(data = boston_redlining_fort, aes(x = long, y = lat, group = group, fill = id2), alpha = 0.8) +
  geom_polygon(data = boston_borders_fort, aes(x= long, y = lat, group = group), fill = NA, color = "gray", alpha = 0.8, size = 0.3) +
  geom_text_repel(data = coordinates, aes(label = neighborhood, x = long, y=lat), color = "black", size = 2, fontface = "bold", point.padding=NA) +
  scale_fill_manual(values = c("#6fd62f", "#0975d8", "#ead723", "#d01313")) +
  labs(title = "Historically Redlined Areas and Today's Boston Neighborhoods", subtitle = "Consequences of 1930s racist federal housing policy affect formally D-graded communities", fill = "HOLC Grade", caption = "Sources: Mapping Inequality & Analyze Boston") +
  theme_void() +
  theme(plot.background = element_rect(fill=NA, color = NA),
        panel.background = element_rect(fill = NA, color = NA),
        plot.title = element_text(size=13, face = "bold", hjust=0.5),
        plot.caption= element_text(size=10),
        plot.subtitle = element_text(size=10, hjust=0.5),
        legend.text = element_text(size = 7),
        legend.title = element_text(size =9)) +
  coord_map()
plot.redlining

#ggsave("img/redlining_map.png", width=8.5, height= 5.5, units="in")

Choropleth 2 & 3: Distribution of POC and low-income populations, gradient for U.S. Census tracts, outline of neighborhood borders

Here is a comparison of the distribution of people of color and low-income populations in Boston.

These maps can be compared to the redlining map. For historically redlined areas and even yellowlined areas, these neighborhoods include more people of color and low-income.

# Choropleth 2: Distribution of POC populations, gradient for U.S. Census tracts, outline of neighborhood borders

#view(climate_ready_sel)
#view(climate_ready_geo_fort)
# Data based on different U.S. Census Tracts
data.climate_ready <- climate_ready_sel %>% 
  select(GEOID10, low_income, POC, tract_population, neighborhood) %>%
  glimpse()
## Rows: 180
## Columns: 5
## $ GEOID10          <dbl> 25025010405, 25025010404, 25025010801, 25025010702, 2…
## $ low_income       <dbl> 1191, 2387, 72, 187, 895, 309, 1024, 1740, 1967, 1560…
## $ POC              <dbl> 1755, 1749, 447, 320, 1039, 517, 1664, 1083, 1123, 19…
## $ tract_population <dbl> 5522, 5817, 2783, 2400, 3173, 3059, 4804, 2791, 4985,…
## $ neighborhood     <chr> "Mission Hill", "Fenway", "Back Bay", "Back Bay", "Fe…
class(data.climate_ready$GEOID10)
## [1] "numeric"
class(climate_ready_geo_fort$id)
## [1] "character"
# convert GEOID10 column to character
data.climate_ready$GEOID10 <- as.character(data.climate_ready$GEOID10)
class(data.climate_ready$GEOID10)
## [1] "character"
# Merge climate_ready datasets (fortified geoJSON: id and csv:GEOID10)
climate_ready_merged <- left_join(climate_ready_geo_fort, data.climate_ready, by=c("id"="GEOID10"))
#view(climate_ready_merged)

# choropleth of POC communities
plot.POC <- ggplot() +
  geom_polygon(data = climate_ready_merged, aes(fill = POC, x = long, y = lat, group=group), alpha = 0.8) +
  geom_polygon(data = boston_borders_fort, aes(x= long, y = lat, group = group),fill=NA, color = "#6821d8", alpha = 0.5, size = 0.3) +
  geom_text_repel(data = coordinates, aes(label = neighborhood, x = long, y=lat), color = "black", fontface = "bold", size = 2, point.padding = NA) +
  scale_fill_gradient(low = "#cfcbd4", high = "#6821d8") +
  labs(title = "Communities of Color Distribution in Boston Neighborhoods", subtitle = "Based on data from 2008-2012 American Community Survey 5-year Estimates, U.S. Census", caption = "Source: Climate Ready Boston", fill = "POC population") +
  theme_void() +
  theme(plot.background = element_rect(fill=NA, color = NA),
        panel.background = element_rect(fill = NA, color = NA),
        plot.title = element_text(size=13, face = "bold", hjust=0.5),
        plot.caption= element_text(size=10),
        plot.subtitle = element_text(size=9, hjust=0.5),
        legend.text = element_text(size = 8),
        legend.title = element_text(size = 10)) +
  coord_map()
plot.POC

#ggsave("img/POC_map.png", width=8.5, height = 5.5, unit="in")
# Choropleth 3: Distribution of low-income populations, gradient for U.S. Census tracts, outline of neighborhood borders
plot.lowIncome <- ggplot() +
  geom_polygon(data = climate_ready_merged, aes(fill = low_income, x = long, y = lat, group=group) ,alpha = 0.8) +
  geom_polygon(data = boston_borders_fort, aes(x= long, y = lat, group = group), fill = NA, color = "#f97100", alpha = 0.5, size = 0.3) +
  geom_text_repel(data = coordinates, aes(label = neighborhood, x = long, y=lat), color = "black", size = 2, fontface = "bold", point.padding = NA) +
  scale_fill_gradient(low = "#ead6c6", high = "#f97100") +
  labs(title = "Low-income Distribution in Boston Neighborhoods", subtitle = "Based on data from 2008-2012 American Community Survey 5-year Estimates, U.S. Census", caption = "Source: Climate Ready Boston", fill = "low-income population") +
  theme_void() +
  theme(plot.background = element_rect(fill=NA, color = NA),
        panel.background = element_rect(fill = NA, color = NA),
        plot.title = element_text(size=13, face="bold", hjust=0.5),
        plot.caption= element_text(size=10),
        plot.subtitle = element_text(size=10, hjust=0.5),
        legend.text = element_text(size = 7),
        legend.title = element_text(size = 9)) +
  coord_map()
plot.lowIncome

#ggsave("img/lowIncome_map.png", width=8.5, height=5.5, unit = "in")

Choropleth 4: Open space areas in Boston, neighborhood borders

One environmental factor to consider is the amount of open space within each neighborhood. Parks and trees within an urban setting are extremely valuable to provide cooling in surrounding areas. People who may not have air-conditioning or other home-cooling options can go outside to these open areas to escape the burden of extreme heat, in event of a heatwave. Yet some communities may have less open space than others.

I obtained a dataset on open space areas of conservation and reservation in Boston, from Analyze Boston. The map below, when compared to the previous maps, shows a correlation of less open spaces in areas that were historically redlined and have greater populations of color and of low-income today.

# Choropleth 4: Open space areas in Boston, neighborhood borders
plot.openspace <- ggplot() +
  geom_polygon(data = open_space_fort, aes(x= long, y = lat, group=group), fill = "forestgreen", alpha = 0.7) +
  geom_polygon(data = boston_borders_fort, aes(x = long, y = lat, group = group), fill = NA, color = "lightgray", alpha = 0.8, size = 0.5) +
  geom_text_repel(data = coordinates, aes(label = neighborhood, x = long, y=lat), color = "black", size = 2, fontface = "bold", point.padding = NA) +
  labs(title = "Open Space in Boston Neighborhoods", caption = "Source: Analyze Boston", subtitle = "Distribution of open green areas can affect urban heat in different neighborhoods") +
  theme_void() +
  theme(plot.background = element_rect(fill=NA, color = NA),
        panel.background = element_rect(fill = NA, color = NA),
        plot.title = element_text(size=13, face = "bold", hjust =0.5),
        plot.caption= element_text(size=10),
        plot.subtitle = element_text(size=10, hjust=0.5)) +
  coord_map()
plot.openspace

#ggsave("img/openSpace_map.png", width=7.5, height = 5, unit = "in")